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Abstract 



The so-called local density approximation plus the multi-orbital mean-field 
Hubbard model (LDA+U) has been implemented within the all-electron pro- 
jector augmented- wave method (PAW), and then used to compute the in- 
sulating antiferromagnetic ground state of NiO and its optical properties. 
The electronic and optical properties have been investigated as a function of 
the Coulomb repulsion parameter U. We find that the value obtained from 
constrained LDA (U = 8 eV) is not the best possible choice, whereas an inter- 
mediate value (U = 5 eV) reproduces the experimental magnetic moment and 
optical properties satisfactorily. At intermediate U, the nature of the band 
gap is a mixture of charge transfer and Mott-Hubbard type, and becomes 
almost purely of the charge-transfer type at higher values of U. This is due 
to the enhancement of the oxygen 2p states near the top of the valence states 
with increasing U value. 
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I. INTRODUCTION 



For many materials, the density functional theory (DFT)S in the local-spin-density ap- 
proximation (LSDA)i provides a good description of their ground-state properties. However, 
problems arise when the DFT-LSDA approach is applied to materials with ions that contain 
incomplete d- or /-shells, such as transition-metal oxides or heavy fermion systems. For ex- 
ample, most transition-metal oxides are wide-gap antiferromagnetic insulators, and the 
DFT-LSDA predicts them to be either metals (FeO and CoO) or small-gap semiconductors 
(MnO and NiO).0 The failure of the DFT-LSDA can be traced to the mean-field character 
of the Kohn-Sham equations as well as to the poor description of strong correlation effects 
within the homogeneous electron gas. The strong correlation effects are responsible for the 
breakdown of the DFT-LSDA description of the electronic structure of these compounds. 
In order to provide a better description of these effects, the Mott-Hubbard picture has been 
introduced.00 

In the Mott-Hubbard picture of NiO, the d-d Coulomb interaction splits the Ni d sub- 
bands into the so-called lower and upper Hubbard bands. The upper Hubbard band has 
mostly Ni 3d 9 character, while the top of the valence band is of 3d 8 character, leading 
to a Mott-Hubbard gap of d-d type. However, O Is x-ray absorptionll! as well as x-ray 
photoemission and bremsstrahlung isochromat spectroscopies! on Li x Ni 1 _ 2 .0 have shown 
that the additional hole has mainly oxygen character. In contrast to the Mott-Hubbard 
model, the energy band gap caused by the Ni 3d correlations is therefore of the charge- 
transfer type between the occupied oxygen 2v and the Ni 3d empty states. 

On the other hand, localized approaches,tjcy such as the local cluster scheme based on 
the configuration interaction method or the Anderson impurity model, in which transition- 
metal ions are treated like an impurity in an oxygen 2p host, predict a well-defined band 
gap of 5.0 eV. However, the oxygen 2p band dispersion observed in angle-resolved photoe- 
mission spectroscopy^ cannot be described by these methods because the lattice effects are 
neglected. 

Several attempts have been made to include the missing correlation effects in DFT- 
LSDA. The generalized gradient approximation (GGA),@ which takes into account the radial 
and angular gradient corrections, can only open a small band gap.@ The self-interaction 
correction 

(sic)iffl 

eliminates the spurious interaction of an electron with itself from the 
conventional DFT-LSDA method. Compared to LSDA, the band gap and the magnetic 
moments are significantly increased. However, the band gap still is too small, and the 
SIC-LSDA method predicts a larger energy band gap for NiO than for FeO and CoO, in 
contradiction to experiment.! The crystal- field orbital polarization introduced by Normanill 
to determine the magnetism and insulating band gap of transition-metal oxides is promising 
but underestimates both the spin magnetic moment and the band gap. 

Another promising approach for correlated materials is the so-called local density ap- 
proximation (LDA) plus the multi-orbital mean-field Hubbard model (LDA+U)liHEI which 
includes the on-site Coulomb interaction in the LSDA Hamiltonian. After adding the on-site 
Coulomb interaction to the LSDA Hamiltonian, the potential becomes spin- and orbital- 
dependent. Because a larger energy cost is associated with fluctuations of the (i-occupancy, 
the orbital-dependent potential reduces the fluctuations of the (i-occupancy, resulting in a 
better justification of a mean-field approach. LDA+U, although it is a mean-field approach, 
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has the advantage of describing both the chemical bonding and the electron-electron inter- 
action. 

The question regarding the best value for the Coulomb repulsion parameter U is, however, 
still under debate. The U parameter for NiO obtained from a constrained LDA calculation 
is about 7 to 8 eV, and this is the value generally used in LDA+U calculations. A similar 
value of U has been obtained from a constrained LDA calculation for bulk FeJH even though 
a much smaller value had been expected because of the metallic screening in Fe. The 
authors arguedi§ that the higher value could be an artifact due to the poor screening within 
the atomic sphere approximation (ASA), and that within a full-potential calculation a much 
smaller value of less than 4 eV would be expected. In contrast, an unpublished full-potential 
linear muffin-tin orbital method calculation by Alouani and Willsll clearly shows that the 
value of U for bulk Fe is even slightly larger than the ASA value. Furthermore, Bulut, 
Scalapino, and WhiteH showed that the renormalization of the Coulomb interaction depends 
on the type of model used. As LDA is not a diagrammatic method, it is not known which 
type of renormalization is the most appropriate for the LDA+U model. 

In this work we shed light on this problem by treating the Coulomb repulsion parameter 
U as adjustable parameter, and by investigating how the electronic and optical properties 
depend on its value. We show that for an intermediate value of U — 5 eV, good agreement 
with the measured ground-state antiferromagnetic magnetic moment and optical properties 
is obtained. We also show that the O 2p character near the top of the valence states is 
enhanced for a larger value of U. Our calculation seems to indicate that the nature of the 
band gap at intermediate U is a mixture of charge transfer and Mott-Hubbard type, and 
that it becomes almost purely of the charge-transfer type for higher values of U. 

Our calculations are based on the projector augmented wave (PAW) method)^ an effi- 
cient all-electron method without shape approximations on the potential or electron density 
to avoid uncertainties due to the ASA approach. Based on a Car-Parrinello-like formalism,H 
the PAW method allows complex relaxations and dynamical properties in strongly corre- 
lated systems to be studied. Our implementation of LDA+U within the PAW method is 
described in detail. Furthermore we discuss possible extensions of the existing method that 
will enhance its applicability. 

The paper is organized as follows: In Section |I| we present those aspects of the PAW 
formalism that are needed for the implementation of the LDA+U method. In Section [ITT] we 
present and discuss the LSDA and LDA+U ground-state properties of NiO, and in Section 
prV| we study the optical properties of NiO, namely the imaginary part of the dielectric 
function, and compare the results to experiment. 



II. FORMALISM 
A. PAW method 

The PAW method developed by one of us0 combines ideas of the pseudo-potential (PP) 
and the linear augmented plane wave (LAPW) methods. It is applicable to all elements of 
the periodic table. The nodal behavior of the wave function is correctly described and, as 
in the PP method, the forces on the ions are easily expressed. 
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In the PAW method, the all-electron (AE) crystal wave function is constructed from a 
pseudo (PS) wave function and atom-like functions localized near the nuclei. The PS wave 
function \^f) coincides with the crystal AE wave function in the interstitial region, i.e. 
outside the atomic regions. Inside the atomic regions Q t , called augmentation regions, the 
wave function is almost atom-like because the effect of the surrounding crystal is small. 
Therefore, a natural choice is to use solutions |0a) of Schrodinger's equation for the isolated 
atom, the so-called AE partial waves, as a basis set for the augmentation region. Here 
A = {t, a,£, m} is a global index for the atom t, the angular momentum £, the magnetic 
quantum number m, and the index a, the energy for which Schrodinger's equation is solved. 

To link the expansion in atom-like functions near the nuclei to the PS wave function, we 
introduce a set of auxiliary functions \4>a), so-called PS partial waves, which are centered 
on the atom and coincide per construction with the corresponding AE partial waves \4>a) 
outside their augmentation regions: 

A (r)=0 A (r) forr^fi t . (1) 

The coefficients ca of the expansions in AE and PS partial waves are chosen such that 
the PS partial wave expansion J2a |0a) c a cancels out the PS wave function inside the 
augmentation region. For this purpose we introduce so-called projector functions (p\ such 
that 

EI#a)(pa| = 1, (2) 

A 

and therefore 

I*) = EI^a>(pa|*> (3) 

A 

for the Hilbert space spanned by the PS partial waves \4>a)- Thus we identify the expansion 
coefficients with ca = (Pa\^)- Equation (0) results in the biorthogonality condition 

(Pa|0a'> = $A,A' ( 4 ) 

for the projector functions, which moreover are chosen to be localized within the correspond- 
ing augmentation region. 

With these conditions, the AE Bloch wave function \l/(r) can be obtained from the PS 
wave function ^(r) as 

*(r) = *(r) + 5>A(r) - ^(r)](p A |*) • (5) 

A 

The first term represents the PS wave function defined over the entire space, which is equal 
to the AE wave function in the interstitial region, and which is expanded in plane waves. The 
second term is the AE partial wave expansion, which describes the correct nodal behavior 
of the wave function in the augmentation region fi^ (r < r*). The third term eliminates the 
spurious contribution of the PS wave function in the augmentation region. 

Note that Eq. (^) holds only approximately if the set of PS partial waves is not entirely 
within the augmentation regions. However, this has the advantage that only those contri- 
butions of the PS wave function will be removed that are also replaced by AE partial waves. 
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As a result, the AE wave function converges rapidly with the number of partial waves used 
and, moreover, it is continuous and different iable for every truncation of the partial-wave 
expansion. 

Expectation values of any sufficiently local operator A are obtained as 

(*L4|*) = + E Wa)((</>aWa') - (4> A \A\4> A ,))(p AI \y) . (6) 

A, A' 

Note that the double sum is diagonal in the site indices t, t' . This equation is exact for a 
complete set of PS partial waves and rapidly attains the converged result if incomplete. The 
PAW method provides the freedom to represent a zero operator in the form 

= - E Wa)(0a|£|<M(pa#> (7) 

A,A' 

by any operator B entirely localized within the augmentation region. Equation (|7|) has the 
same range of validity as Eq. (|5]) does and allows a further acceleration of the convergence by 
using a well-chosen operator B and adding the corresponding zero operator to the expression 
for the expectation value. 



B. The LDA+U total energy functional 

For transition-metal oxides, the <i-orbitals are well localized and keep a strong atom- 
like character. Even though LDA provides a good approximation for the average Coulomb 
energy of the d-d interactions, it fails to describe correctly the strong Coulomb and ex- 
change interaction between electrons in the same d-shell. The main intention of LDA+U 
is to identify these atomic orbitals and to describe their electronic interactions as strongly 
correlated states. The other orbitals are delocalized and considered to be properly described 
by the LDA. The procedure is to eliminate the averaged LDA energy contribution of these 
atom-like orbitals from the LDA total energy functional E LDA , and to add an orbital- and 
spin-dependent correction. The total energy within the LDA+U method then has the form 

e = £ LDA + \ E E (xh xilKelx'; x'HHT 

t,u i,j,k,l 



2 

t,a i,j,k,l 



E 



^ _ _/ ^ ST 



(7.(7 



(8) 



where N t,cr = YU n i1 i s the average occupation of the <i-shell for each spin direction as 



obtained from the (i-orbital occupancies n*|J . U and J are the Coulomb self-energy and 
the exchange parameter, respectively. The expressions Xk\Vee\Xj] x\) are the four-center 
matrix elements of the screened Coulomb interaction V ee . An additional requirement of the 
LDA+U approach is that the additional energy is applied only to the valence electrons, 
which are re-optimized while constrained to remain orthogonal to the core states. 
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C. Orbital occupations 



The rationale behind the definition of the orbital occupations is to project the density 
matrix onto one consisting of a restricted set of localized orbitals where each orbital 
index i stands for an angular momentum and t for an atomic site. Here we focus in particular 
on the o?-orbitals of Ni. 

The projected density matrix operator has the form 

^ = £lx'KJ<4l- ( 9 ) 

The density matrix elements n\ ^ in the localized basis set can be obtained from the condition 
that the expectation values of a set of projection operators P s are optimally reproduced by 
the localized basis set, i.e. the weighted square deviation 



E 



E(*nl^l*n ,<7 )-E^;J(x*l^lx5) 

n,k i,j 



of the projections should be minimized. 

Here we choose the projection operators acting on the Ni <i-orbitals in analogy to previous 
implementations^ as 

P^ >m ,(r, r') = e nt (r)5(\r' - R t | - |r - R*|)^ Jr)y d * m ,(r') , (10) 

where the site index t refers to a particular Ni site, and \Y dm ) is the spherical harmonic cen- 
tered at site t. The step functions 9n t (r) are unity for \r — R t \ < r l c and zero otherwise. Y djTn 
are the spherical harmonics for i = 2. (Note that (x\P\x') = I dr J d r 'x{ r )P{ r , r ')x'{ r ')-) 
The radius r* for Ni corresponds to an atomic sphere radius, and we have chosen r* = 2.1 
ao- The localized orbitals Ixln) have been chosen to be identical to those of the spheridized, 
non-spin-polarized atoms. We used an atom in the 3d 8 As 2 configuration of Ni. For such a 
localized Ni <i-orbital Ixln) with magnetic quantum number m, we obtain 

I (Xml-^TT^mlXm) (Xm I ^fit (?) I Xm) ' (H) 

independent of m and t. Thus we obtain the orbital occupations 

1 

n,k 



t,a 
m.m' 



jY.^Pl^Y) (12) 

' 1 . 



from the minimal condition for F, which in this case are independent of the weights w s . 

In the PAW method we obtain the orbital occupations directly from the PS wave func- 
tions \^n' a ) aS 



1 

m.m' 



7 E<^'1^,J*!r>, (13) 

k,n 



using the pseudo- version 
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Pl,m' = E \PA)(<f>A\PL, m >\M<PA'\ + A ^m,m' (14) 



A,A' 

of the projection operator P^m'- Consistent with the PAW formalism, the small correction 

A^U = Pl,m> - E \PA){k\PL,m¥A')(PA'\ (15) 
A,A' 

is ignored in the present calculations because it can be considered the pseudo-version of the 
zero operator. 

Note that the present choice for the projection operator may not be the most elegant 
form. Here we suggest another version using projector functions (qj\. The construction is 
analogous that of the projector functions in the PAW method, but the purpose is different: 
namely, the goal is to decompose a wave function into local orbitals such that 

Elx*>teW = l*> (is) 

i,t 

holds for any wave function that can be represented exactly as a superposition of the 
local orbitals |x*). This requirement is fulfilled whenever (<7*|x*) = Sij5 tt t'- The projector 
functions can be obtained from a set of functions |/*), which are linearly independent and 
localized within an augmentation region as 

E</ii4'><<?il = <//l (17) 
j,t' 

after multiplication with the inverse of (/*|x*-)- If we now use 

Pt j = \ci)(q t j \, (18) 
we obtain the matrix elements in the representation of localized orbitals directly as 



k,n 



It is now a simple matter to obtain the corresponding pseudo-projection operator 

Kj = E |PA><0A|$<gJ|0A'><PA'| • (20) 



This expression has a number of advantages: (i) It is the most general expression because 
any operator can be expressed in a fully separable expression,^! and, unlike Eq. fllPf) , it is 
not limited to a semi-local form, (ii) The step function in Eq. (pi]) introduces a potential 
that is discontinuous, whereas the operator suggested here can be constructed such that the 
resulting potential is continuous and differentiable, while still being fully localized in a well- 
defined region, (iii) The separable form of the operator renders the evaluation numerically 
convenient. In particular, it is conceivable to evaluate AP explicitly if necessary, whereas the 
corresponding operation using the semi-local projection operator would be computationally 
prohibitive, (iv) We can systematically construct projection operators that specifically act 
on a well-defined portion of the Hilbert space: By enlarging the set of functions to 
include also orbitals that should not be represented in the restricted set, they are excluded 
by simply setting the corresponding matrix elements of the density matrix Nij to zero. 
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D. Coulomb and exchange parameters 



The four-center integrals used in the expression of the LDA+U total energy are defined 

as 

(xhxjlVeelXklxl) = J d ^J dr.xT^xfi^v^r^xK^xK^) , (21) 

where v ee (r, r') is the screened Coulomb interaction between two electrons. 

If we choose localized <i-orbitals that are described by an atomic <i-wave function Xm( r ) = 
Xd(\ r — Rt\)Yd, m ( r — R-t) with magnetic quantum number m (where Y^ m are the spherical 
harmonics for t = 2) and furthermore assume that the static dielectric function e is constant 
in space, we can exploit the multipole expansion of l/|ri — r2 1: 

Mri, r 2 ) = = - ± ^hii E YtJrxYLfr) ■ (22) 

e|ri-r 2 | e e=Q M + 1 r> m= _ e 

Here r< and r> denote the smallest and largest values of r\ and r 2 , respectively Under these 
assumptions we can transform Eq. ( pij ) into 

OO A +£ 

(xi;X3lKe|X2;xi> = £^-7 E (ii,m 1 \Y i} j£ 2 ,m 2 )(£ 3 ,m 3 \Y e *j£ 4 ,m i )F i , (23) 

£=0 m=—£ 

where (£i,rrii) are the angular momenta quantum numbers of (£, m\Y£" jm "\£, ml) the 
Gaunt coefficients, and F the so-called screened Slater's integrals. Because of the special 
properties of the Gaunt coefficients, only F°,F 2 , and F A contribute to the Coulomb integrals: 

T /*oo /*oo rp£ 

F l = ~ drj dr 2 rfrlxKr^xKr,)^ . (24) 
e Jo Jo r> 

The parameters f/ and J are identified with averages of the Coulomb and exchange 
integrals, which are related to the Slater integrals F°,F 2 , and F 4 by the properties of 
Clebsch-Gordan coefficients, Eq. (|23|): 

1 t t t t o 

^ (2£ + l) 2 -^m' I ^ee | Xm ) Xm') F 5 (25) 



iii.rn 



J= , n E (xL;x^lK e |x^;xL) = (^ 2 + ^ 4 )/i4, (26) 



(2£)(2£- 



The dielectric constant and therefore the Coulomb and exchange parameters {7 and J are not 
known a priori. Usually, they are obtained from a constrained DFT calculation.il Here, we 
are interested in how sensitively the results depend on the choice of Coulomb parameters, and 
which Coulomb parameters will provide the best agreement with reality as probed by optical 
absorption. Therefore, we adopt the general form for the four-center integrals as function 
of the U and J suggested by the arguments provided above, and perform calculations for 
different U values, namely U = 5 eVx/ 2 and U = 8 eVxI 2 . In the following we suppose 
that I 2 = 1. Because the results are fairly insensitive to the exchange parameter J, we 
have adopted J = 0.95 eVx/ 2 from previous constrained LDA calculations.il The third 
mandatory relation is obtained from the work of DeGroot et al.,0 who determined the ratio 
F A jF 2 for transition-metal oxides to be between 0.62 and 0.63. We therefore adopt a ratio 
F 4 /F 2 = 0.625. 



S 



E. Hamiltonian 



The pseudo-Hamiltonian operator 

H a = H L a DA + H", (27) 

which acts on the PS wave functions, is obtained as the derivative of the total en- 
ergy functional with respect to the two-center pseudo-density matrix operator p = 
En,k \ ^n' a )fn' a (^n a \- Tne non-LDA contribution of the LDA+U Hamiltonian is then ob- 
tained as the product of the derivative of the non-LDA contribution to the total energy 
and the projection operator P, which is the derivative of the occupation with respect to the 
two-center pseudo-density matrix operator. Thus we obtain 

t,m,m' 

where 



v: 



t,cr 

mi ,m2 



rri3,iri4 



= V" (y l * y* W Iy* -y* )n t ' 

/ , \Ami' Am.3 1 eelAm2' Am,4/ m 
m3,rrt4 



7(13,7714 



n 



E 



*7(iV 



(29) 



The LDA contribution of the pseudo-Hamiltonian has the usual formS 



H. 



LDA 



V 



— + v+ J2 \Pm) 



Ai,A 2 



V 2 



V 2 ~ 



(PA 2 | • (30) 



For the LDA+U calculation, we first performed a self-consistent LSDA calculation using 
the all-electron PAW method, and then used the self-consistent potential to construct the 
LSDA Hamiltonian for a large number of k-points in the Brillouin zone. Next, the Hubbard 
correction is added to the LSDA Hamiltonian as given by Eq. ([27]), and the new Hamiltonian 
iterated until the occupation numbers n^ m , have converged. We did not use the so-called 
second-variation procedure for self-consistent LDA+Uj^l in which the LSDA potential is 
updated, because it failed to yield any improvement on calculations similar to ours.il 



III. GROUND STATE OF NiO 

A. LSDA ground state 

The ground state of NiO has been calculated using the PAW method, and the density 
of states (DOS) is calculated from the self-consistent PAW potential using the tetrahedron 
method for the Brillouin-zone integration. eil Figure [I] presents the atom-resolved DOS in 
the augmentation region. As can be seen, LSDA produced an antiferromagnetic insulating 
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ground state with a small band gap. The LSDA magnetic moment is 0.95 and mainly 
due to the Ni £g band splitting. This value is much smaller than the experimental value 
(1.64-1.90 /is, Refs. f42| , f43|) . The calculated band gap of about 0.1 eV is also much smaller 
than the experimental one (3.0-4.0 eV). The most interesting feature of our LSDA DOS 
is that the ci-states of Ni dominate the region in the vicinity of the band gap, and that 
the top of the valence state is of Ni eg type for the first and Ni i29 for the second spin. 
This electronic structure suggests that the band gap is of Mott-Hubbard type. Hence, this 
LSDA picture of NiO disagrees completely with experiment. It is surprising that the quasi- 
particle calculation within the so-called GW approximation, performed by Aryasetiawan 
and Gunnarsson,El produced results qualitatively similar to LSDA except for an increased 
band gap of 6 eV and an increased magnetic moment of 1.6 /ie. The quantitative change 
is the reduction of the O 2p bandwidth bvalmost 1 eV. However, a recent self-consistent 
model GW calculation by Massida et a/.,c3 in which the dynamic effects were neglected, 
produced other results than those of Aryasetiawan and GunnarssonS It was argued by 
Massida al. that the results of the former GW calculation are not quite self-consistent, 
presumably because of the additional non-local ad- hoc potential that is adjusted to the GW 
calculation in each self-consistent step. The main difference between the two reported GW 
model calculations is that the latter calculation^ produced (i) a spreading of the Ni <i-states 
over the entire valence bandwidth, (ii) a vanishing gap between the O 2p and Ni 3d and, 
most importantly, (iii) an enhancement of O 2p states at the valence band maximum. The 
latter effect attributes the origin of the band gap mainly to a charge transfer gap because 
this gap is now between the O 2p and Ni 3d conduction states. Next, we will show that this 
latter finding is in agreement with the results of the LDA+U model. 



B. LDA+U ground state 

We have used our implementation of the LDA+U model to determine the ground-state 
electronic structure of NiO. Although it is common practice to use the U extracted from 
a constrained LDA calculation, we adopt a different point of view here. As stated in the 
introduction and in ag reement with recent results reported in the literature^ believe that 
the value of U extracted from constrained LDA is not the best possible choice. Therefore 
we have determined the electronic structure of NiO for an intermediate U of 5 eV as well as 
for a larger value of 8 eV. 

Figure ^ shows our calculated LDA+U DOS for U = 5 and 8 eV. The energy band gap 
is found to be 2.8 and 4.1 eV, respectively. The total antiferromagnetic spin moment is 
1.73 and 1.83 /zg, respectively. Our DOS obtained for U = 8 eV agrees well with previous 
LDA+U calculations.tlSEiElci However, the DOS for U = 5 eV is in better agreement with 
the GW model calculation of Massida et a/.B The top of valence band is reinforced by the 
O 2p states, rendering the band gap a mixture of charge-transfer type and Ni <i — *> <i-like 
excitations. In agreement with previous LDA+U calculations, the spin majority Ni e g states 
are pushed towards lower energies, and the energy difference between and is about 11 
eV for U = 8 eV and 8.6 eV for U = 5 eV. Here again the GW model yields a value of about 
9 eV for this splitting, in good agreement with our results for U = 5 eV. 
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IV. DIELECTRIC FUNCTION 



For insulators, the imaginary part of the macroscopic dielectric function is obtained 
within the random phase approximation (RPA) in the long wave length limit without local- 
field effects as@ 

e 2 (u) = hm V E S l M 2-J 2 /«,k(l - - % - M • (31) 

q v,c k " 

Here M^' k _ q are the interband transition matrix elements, /^k is the zero-temperature Fermi 
distribution, Q is the cell volume, c denotes the conduction-band and v the valence-band 
index. In the case of a local potential, the interband transition matrix elements are given 
by 

r c,k _ Q /vTvki |,T,k 



limM^_ q = 2 <| p |v^ , (32) 



where the matrix elements (\l/ k |p|\l/ k ) are calculated using the PAW crystal wave functions 
\l/ k described by Eq. (|): 

(* k |p|* k ) = (^ k | P |^ k ) + E<*«Ipa) [(0a|p|0A') - (0a|p|<M] (pa'|*J) , (33) 

A, A' 

In the most general case, i.e. where the potential is nonlocal as in LDA+U, a nonlocal con- 
tribution has to be added to the interband transition matrix elements.^ The full derivation 
is given in the Appendix by Eq. ( |A~6f) . For NiO the nonlocal contribution to the matrix 
elements is found to be small, i.e. of a few percent. 

Figure ^| shows the imaginary part of the dielectric function calculated within the LSDA. 
The resulting optical spectrum is not in agreement with experiment (see Fig. ^). In particu- 
lar, the optical gap is considerably underestimated and the first structure has a much higher 
intensity compared to experiment. Conversely, for U = 5 eV, our calculated imaginary part 
of the dielectric function within the LDA+U is in a better agreement with experiment, as 
shown in Fig. |]. The optical band gap and the oscillator strength of the first excitation 
peak are in excellent agreement with experiment. However, at higher photon energies the 
agreement with experiment is only qualitative, which is not expected owing to the mean- 
field approximation of this simple model. A much higher value of U, i.e. 8 eV, produces a 
much larger optical gap in contrast to experiment. In agreement with our conclusion that 
a much smaller value of U is required to describe NiO, Dudarev and coworkers^ also found 
that U = 6.2 eV reproduces the lattice parameter and the measured electron-energy loss 
spectra. It is too early to draw a definitive conclusion about the excited states of NiO, as 
our LDA+U model is a mean-field-like model in which excitonic effects are not included. To 
our knowledge, no calculations of excitonic effects have so far been attempted. Our results 
represent the first investigation of the low-lying excited states of NiO that considers all the 
subtleties of chemical bonding and strong electron-electron interaction. 

The interband transitions are responsible for the first structure in the optical spectrum 
of NiO, located between 4.1 and 4.4 eV. We found that 40.2% of the contribution results 
from the transition from band 15 (second highest occupied band) to band 17 (lowest empty 
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band), 36.2% from the transition 16— >18, and 15.9% from the interband transition 16 — >17. 
To analyze the character of the initial and final states of the interband transitions, the band 
structure dispersion along some of the high-symmetry directions is shown in Fig. [5| together 
with the DOS of the states that give rise to the first optical peak. The arrows between the 
parallel bands indicate the interband transitions from the initial to the final state responsible 
for the first structure in the optical spectrum. Figure || shows the charge density plot of 
the initial and final states \E^ ,<T for bands 16 (highest occupied band) and 18 (second lowest 
empty band) at point k = (jf^? Mo"' o l° ca t e d between the high-symmetry points K and 
U, where the optical matrix element value is among the largest. It also shows that the 
initial state is of mixed O 2p and Ni 3c? character, whereas the final state mainly is of Ni e g 
character. The optical transitions then are between the O 2p and the Ni 3d states, resulting 
in excitations of the charge-transfer type. Our interpretation differs from that of Fujimori 
and Minami,@ who used a configuration interaction within the metal-ligand cluster and 
claimed that the d — > d charge-transfer transitions are the origin of fundamental edge. The 
drawback of the cluster calculation is that in reality the oxygen 2p orbitals are delocalized 
and therefore not well described in a small cluster. On the other hand, an earlier, band- 
structure-based interpretation by Messick aud coworkers! assigned the peak to one electron 
interband transitions associated with Ni 3d to the Ni 4s states. This interpretation is not 
correct either because the Ni 4s are far above the top of the valence states (greater than 
6 eV), and only quadrupolar interband transitions are permitted between the 3d and 4s, 
which substantially reduces the peak intensity. 



V. CONCLUSION 

We have presented a new implementation of LDA+U model based on the PAW method,!! 
which is an all-electron method without any shape approximation to the potential or the 
charge density. We tested the method on NiO and obtained results that are in good agree- 
ment with previous LDA+U calculations and a recent GW model calculation. In particular, 
we obtained the correct antiferromagnetic insulating ground state of NiO. 

We discussed the results in terms of the strength of the Hubbard interaction U. The 
optimum value of U has been determined by comparison with the experimental dielectric 
function as well as with the ground state properties. We observed a large enhancement 
of the O 2p character at the top of the valence state, resulting in a more charge-transfer 
than Ni d — > d LSDA -type band gap. The calculated antiferromagnetic moment is in good 
agreement with experiment. 

Our calculated dielectric function for an intermediate value of U, namely 5 eV, is in 
good agreement with experiment. The low- lying, strong structure in the optical spectrum 
has been assigned to an interband transition from O 2p states at top of the valence band to 
the Ni e g states at the conduction-band bottom. Hence the origin of the first optical peak 
is due to a charge-transfer excitation. 

Our calculation is supported by a recent LDA+U calculation by Dudarev and 
coworkers j^l who also argue that a much smaller value of U than the one obtained from 
constrained LDA calculation is needed to describe the electron energy loss spectra and 
the equilibrium lattice parameter. It should be interesting to apply this method to other 
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transition-metal oxides and check the applicability of LSDA+U for producing excitation 
energies. 
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APPENDIX A: OPTICAL TRANSITION MATRIX ELEMENTS IN LDA+U 



The interband transition matrix elements for a given Hamiltonian are obtained as follows: 



Kt q =( VI/ ^ q |e- Jqr |vl/ c , k ) 



(fr„,k- q |(e»,,k-q ~ e c ,k)e~ iqr |^ c , k ) 

^?j,k— q ^c,k 



(^ ik _ q |Fe-^ r - e-^tfl^O _ (^,k-q|[^,e- iqr ]_|^ c , k ) 



&»,k— q ^c,k 



^i),k— q ^c,k 



(Al) 



The commutator involving of the LDA contribution iJ LDA to the Hamiltonian H = H LDA + 
H u is 



[# LDA ,e~ iqr ]- = -^[V 2 ,e- iqr ] 



-2iqV + g 2 ) = -qp + 0(g) 2 . 



(A2) 



The quadratic and higher-order terms in q can be ignored in the long-wavelength limit 
appropriate for optical transitions. The commutator involving the non-LDA part H u is 
obtained as 



Vl 



t,m,m' 



(A3) 



Next we use the relation which holds for the special form of the projector operator presented 
in Eq. (0) 



-^m",m'l0<^,m,a) _ @£lt ( r ) \4>t,l,m" ,a)bt,t'$i,2$m,r 



(A4) 



i.e. the projection operator changes the magnetic quantum number of a specific d-\ike partial 
wave from m to m" and truncates it beyond the atomic sphere Q t - Hence we obtain 



(A5) 



Finally, we obtain the expression of the matrix elements for the dipole transition, with 
the PAW LDA+U formalism: 



lim M'„ 



c,k 

i>,k— q 



q 



{^v,k-q e c,kJ [ A, A' 

V. 



|p|<M - (0a|p|0A') 



+ i E -^((^V^Al^Wr^A') - (0A|^(r)r|P4 /m ,0 A ')) 



(pA'l^ck) • (A6) 



The difference between wave functions l^.k) and l^k-q) has been ignored because it only 
contributes to terms that are proportional to q 2 , which are ignored in the long- wavelength 
limit. 
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TABLES 



TABLE I. Magnetic moment and band gap of NiO within LDA, LDA+U (U = 5 eV) and 
LDA+U ({7 = 8 eV). The radii of the augmentation regions are 1.7 and 2.1 a.u. for the oxygen 
and the nickel atoms, respectively. 





LDA 


LDA+U 


LDA+U 


Expt. 






U = 5 eV 


{7 = 8 eV 




m (fi B ) 


0.95 


1.73 


1.83 


1.64 (Ref. H) - 1.9 (Ref. ||| 


gap (eV) 


0.1 


2.8 


4.1 


3 - 4.4 
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FIGURES 




FIG. 1. Atom-resolved antiferromagnetic density of states (DOS, in states per unit cell per eV) 
of NiO calculated with LSDA. The band gap is, with about 0.1 eV, significantly underestimated. 
The spin magnetic moment is 0.95 \ib- 
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FIG. 2. Atom-resolved density of states (DOS, in states per unit cell per eV) of NiO calculated 
with LDA+U for (a) U = 5 eV and (b) U = 8 eV . The calculated band gap is significantly improved 
and is ~ 2.8 eV (U = 5 eV) and 4.1 eV (U = 8 eV). The antiferromagnetic total magnetic moment 
of 1.73 (is f° r {U = 5 eV) and 1.83 /i£ for (U = 8 eV) is also in good agreement with experiment. 
The top of the valence state is now O 2p-like, thus producing a mixed charge-transfer-type band 
gap. 
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FIG. 3. LSDA calculated imaginary part of the dielectric function €2(0;) of NiO. Compared to 
experiment, the optical energy gap is underestimated and the first excitation peak is very intense. 
This intense peak is due to interband transitions from the top of the valence band of O 2p character 
to the bottom of the conduction band of type Ni e g character. 
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FIG. 4. Calculated LDA+U imaginary part of the dielectric function of NiO convoluted with a 
Gaussian lifetime at half maximum of 0.3 eV for U = 5 eV (solid curve) and for [7 = 8 eV (dashed 
curve) compared to the experimental optical spectrum! (dashed-dotted curve). The agreement 
with experiment is much better for U = 5 eV. Most of the interband transitions giving rise to the 
first peak are from the top of the valence band of O 2 character to the bottom of the conduction 
band of Ni e g character. At higher photon energy, the disagreement becomes stronger. 
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FIG. 5. Left panel: band structure of NiO calculated with LDA+U for U = 5 eV along the 
high-symmetry directions TZ, ZK, and KU. Right panel: density of states (DOS, in states per unit 
cell per eV) of the initial and final states responsible for the first peak of the imaginary part of the 
dielectric function. The interband transitions between parallel bands giving rise to the first optical 
peak are indicated by arrows. 
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FIG. 6. Initial (top) and final (bottom) charge densities of one spin direction at the point 

k = (iM' Ms' ^W^)a located between the high-symmetry points K and U of the topmost valence 
band (16) and the second lowest conduction band (18), where the optical matrix element value is 
among the largest. The interband transitions giving rise to the first peak of the optical spectrum 
of NiO are due to allowed electronic transition between initial O 2p and Ni e g final states as shown 
in the plot. 
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